<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>shseqr.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SHSEQR.1"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
     $                   LDZ, WORK, LWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
      CHARACTER          COMPZ, JOB
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               H( LDH, * ), WI( * ), WORK( * ), WR( * ),
     $                   Z( LDZ, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     Purpose
</span><span class="comment">*</span><span class="comment">     =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     <a name="SHSEQR.19"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a> computes the eigenvalues of a Hessenberg matrix H
</span><span class="comment">*</span><span class="comment">     and, optionally, the matrices T and Z from the Schur decomposition
</span><span class="comment">*</span><span class="comment">     H = Z T Z**T, where T is an upper quasi-triangular matrix (the
</span><span class="comment">*</span><span class="comment">     Schur form), and Z is the orthogonal matrix of Schur vectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Optionally Z may be postmultiplied into an input orthogonal
</span><span class="comment">*</span><span class="comment">     matrix Q so that this routine can give the Schur factorization
</span><span class="comment">*</span><span class="comment">     of a matrix A which has been reduced to the Hessenberg form H
</span><span class="comment">*</span><span class="comment">     by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Arguments
</span><span class="comment">*</span><span class="comment">     =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     JOB   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">           = 'E':  compute eigenvalues only;
</span><span class="comment">*</span><span class="comment">           = 'S':  compute eigenvalues and the Schur form T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     COMPZ (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">           = 'N':  no Schur vectors are computed;
</span><span class="comment">*</span><span class="comment">           = 'I':  Z is initialized to the unit matrix and the matrix Z
</span><span class="comment">*</span><span class="comment">                   of Schur vectors of H is returned;
</span><span class="comment">*</span><span class="comment">           = 'V':  Z must contain an orthogonal matrix Q on entry, and
</span><span class="comment">*</span><span class="comment">                   the product Q*Z is returned.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     N     (input) INTEGER
</span><span class="comment">*</span><span class="comment">           The order of the matrix H.  N .GE. 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ILO   (input) INTEGER
</span><span class="comment">*</span><span class="comment">     IHI   (input) INTEGER
</span><span class="comment">*</span><span class="comment">           It is assumed that H is already upper triangular in rows
</span><span class="comment">*</span><span class="comment">           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
</span><span class="comment">*</span><span class="comment">           set by a previous call to <a name="SGEBAL.50"></a><a href="sgebal.f.html#SGEBAL.1">SGEBAL</a>, and then passed to <a name="SGEHRD.50"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>
</span><span class="comment">*</span><span class="comment">           when the matrix output by <a name="SGEBAL.51"></a><a href="sgebal.f.html#SGEBAL.1">SGEBAL</a> is reduced to Hessenberg
</span><span class="comment">*</span><span class="comment">           form. Otherwise ILO and IHI should be set to 1 and N
</span><span class="comment">*</span><span class="comment">           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
</span><span class="comment">*</span><span class="comment">           If N = 0, then ILO = 1 and IHI = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     H     (input/output) REAL array, dimension (LDH,N)
</span><span class="comment">*</span><span class="comment">           On entry, the upper Hessenberg matrix H.
</span><span class="comment">*</span><span class="comment">           On exit, if INFO = 0 and JOB = 'S', then H contains the
</span><span class="comment">*</span><span class="comment">           upper quasi-triangular matrix T from the Schur decomposition
</span><span class="comment">*</span><span class="comment">           (the Schur form); 2-by-2 diagonal blocks (corresponding to
</span><span class="comment">*</span><span class="comment">           complex conjugate pairs of eigenvalues) are returned in
</span><span class="comment">*</span><span class="comment">           standard form, with H(i,i) = H(i+1,i+1) and
</span><span class="comment">*</span><span class="comment">           H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
</span><span class="comment">*</span><span class="comment">           contents of H are unspecified on exit.  (The output value of
</span><span class="comment">*</span><span class="comment">           H when INFO.GT.0 is given under the description of INFO
</span><span class="comment">*</span><span class="comment">           below.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Unlike earlier versions of <a name="SHSEQR.68"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>, this subroutine may
</span><span class="comment">*</span><span class="comment">           explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
</span><span class="comment">*</span><span class="comment">           or j = IHI+1, IHI+2, ... N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     LDH   (input) INTEGER
</span><span class="comment">*</span><span class="comment">           The leading dimension of the array H. LDH .GE. max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     WR    (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">     WI    (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">           The real and imaginary parts, respectively, of the computed
</span><span class="comment">*</span><span class="comment">           eigenvalues. If two eigenvalues are computed as a complex
</span><span class="comment">*</span><span class="comment">           conjugate pair, they are stored in consecutive elements of
</span><span class="comment">*</span><span class="comment">           WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
</span><span class="comment">*</span><span class="comment">           WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
</span><span class="comment">*</span><span class="comment">           the same order as on the diagonal of the Schur form returned
</span><span class="comment">*</span><span class="comment">           in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
</span><span class="comment">*</span><span class="comment">           diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
</span><span class="comment">*</span><span class="comment">           WI(i+1) = -WI(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Z     (input/output) REAL array, dimension (LDZ,N)
</span><span class="comment">*</span><span class="comment">           If COMPZ = 'N', Z is not referenced.
</span><span class="comment">*</span><span class="comment">           If COMPZ = 'I', on entry Z need not be set and on exit,
</span><span class="comment">*</span><span class="comment">           if INFO = 0, Z contains the orthogonal matrix Z of the Schur
</span><span class="comment">*</span><span class="comment">           vectors of H.  If COMPZ = 'V', on entry Z must contain an
</span><span class="comment">*</span><span class="comment">           N-by-N matrix Q, which is assumed to be equal to the unit
</span><span class="comment">*</span><span class="comment">           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
</span><span class="comment">*</span><span class="comment">           if INFO = 0, Z contains Q*Z.
</span><span class="comment">*</span><span class="comment">           Normally Q is the orthogonal matrix generated by <a name="SORGHR.95"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>
</span><span class="comment">*</span><span class="comment">           after the call to <a name="SGEHRD.96"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a> which formed the Hessenberg matrix
</span><span class="comment">*</span><span class="comment">           H. (The output value of Z when INFO.GT.0 is given under
</span><span class="comment">*</span><span class="comment">           the description of INFO below.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     LDZ   (input) INTEGER
</span><span class="comment">*</span><span class="comment">           The leading dimension of the array Z.  if COMPZ = 'I' or
</span><span class="comment">*</span><span class="comment">           COMPZ = 'V', then LDZ.GE.MAX(1,N).  Otherwize, LDZ.GE.1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     WORK  (workspace/output) REAL array, dimension (LWORK)
</span><span class="comment">*</span><span class="comment">           On exit, if INFO = 0, WORK(1) returns an estimate of
</span><span class="comment">*</span><span class="comment">           the optimal value for LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     LWORK (input) INTEGER
</span><span class="comment">*</span><span class="comment">           The dimension of the array WORK.  LWORK .GE. max(1,N)
</span><span class="comment">*</span><span class="comment">           is sufficient, but LWORK typically as large as 6*N may
</span><span class="comment">*</span><span class="comment">           be required for optimal performance.  A workspace query
</span><span class="comment">*</span><span class="comment">           to determine the optimal workspace size is recommended.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If LWORK = -1, then <a name="SHSEQR.114"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a> does a workspace query.
</span><span class="comment">*</span><span class="comment">           In this case, <a name="SHSEQR.115"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a> checks the input parameters and
</span><span class="comment">*</span><span class="comment">           estimates the optimal workspace size for the given
</span><span class="comment">*</span><span class="comment">           values of N, ILO and IHI.  The estimate is returned
</span><span class="comment">*</span><span class="comment">           in WORK(1).  No error message related to LWORK is
</span><span class="comment">*</span><span class="comment">           issued by <a name="XERBLA.119"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.  Neither H nor Z are accessed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     INFO  (output) INTEGER
</span><span class="comment">*</span><span class="comment">             =  0:  successful exit
</span><span class="comment">*</span><span class="comment">           .LT. 0:  if INFO = -i, the i-th argument had an illegal
</span><span class="comment">*</span><span class="comment">                    value
</span><span class="comment">*</span><span class="comment">           .GT. 0:  if INFO = i, <a name="SHSEQR.126"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a> failed to compute all of
</span><span class="comment">*</span><span class="comment">                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
</span><span class="comment">*</span><span class="comment">                and WI contain those eigenvalues which have been
</span><span class="comment">*</span><span class="comment">                successfully computed.  (Failures are rare.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                If INFO .GT. 0 and JOB = 'E', then on exit, the
</span><span class="comment">*</span><span class="comment">                remaining unconverged eigenvalues are the eigen-
</span><span class="comment">*</span><span class="comment">                values of the upper Hessenberg matrix rows and
</span><span class="comment">*</span><span class="comment">                columns ILO through INFO of the final, output
</span><span class="comment">*</span><span class="comment">                value of H.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                If INFO .GT. 0 and JOB   = 'S', then on exit
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           (*)  (initial value of H)*U  = U*(final value of H)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                where U is an orthogonal matrix.  The final
</span><span class="comment">*</span><span class="comment">                value of H is upper Hessenberg and quasi-triangular
</span><span class="comment">*</span><span class="comment">                in rows and columns INFO+1 through IHI.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                If INFO .GT. 0 and COMPZ = 'V', then on exit
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                  (final value of Z)  =  (initial value of Z)*U
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                where U is the orthogonal matrix in (*) (regard-
</span><span class="comment">*</span><span class="comment">                less of the value of JOB.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                If INFO .GT. 0 and COMPZ = 'I', then on exit
</span><span class="comment">*</span><span class="comment">                      (final value of Z)  = U
</span><span class="comment">*</span><span class="comment">                where U is the orthogonal matrix in (*) (regard-
</span><span class="comment">*</span><span class="comment">                less of the value of JOB.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                If INFO .GT. 0 and COMPZ = 'N', then Z is not
</span><span class="comment">*</span><span class="comment">                accessed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ================================================================
</span><span class="comment">*</span><span class="comment">             Default values supplied by
</span><span class="comment">*</span><span class="comment">             <a name="ILAENV.162"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>(ISPEC,'<a name="SHSEQR.162"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
</span><span class="comment">*</span><span class="comment">             It is suggested that these defaults be adjusted in order
</span><span class="comment">*</span><span class="comment">             to attain best performance in each particular
</span><span class="comment">*</span><span class="comment">             computational environment.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            ISPEC=1:  The <a name="SLAHQR.167"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> vs <a name="SLAQR0.167"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a> crossover point.
</span><span class="comment">*</span><span class="comment">                      Default: 75. (Must be at least 11.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            ISPEC=2:  Recommended deflation window size.
</span><span class="comment">*</span><span class="comment">                      This depends on ILO, IHI and NS.  NS is the
</span><span class="comment">*</span><span class="comment">                      number of simultaneous shifts returned
</span><span class="comment">*</span><span class="comment">                      by <a name="ILAENV.173"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>(ISPEC=4).  (See ISPEC=4 below.)
</span><span class="comment">*</span><span class="comment">                      The default for (IHI-ILO+1).LE.500 is NS.
</span><span class="comment">*</span><span class="comment">                      The default for (IHI-ILO+1).GT.500 is 3*NS/2.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            ISPEC=3:  Nibble crossover point. (See <a name="ILAENV.177"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> for
</span><span class="comment">*</span><span class="comment">                      details.)  Default: 14% of deflation window
</span><span class="comment">*</span><span class="comment">                      size.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            ISPEC=4:  Number of simultaneous shifts, NS, in
</span><span class="comment">*</span><span class="comment">                      a multi-shift QR iteration.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      If IHI-ILO+1 is ...
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      greater than      ...but less    ... the
</span><span class="comment">*</span><span class="comment">                      or equal to ...      than        default is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                           1               30          NS -   2(+)
</span><span class="comment">*</span><span class="comment">                          30               60          NS -   4(+)
</span><span class="comment">*</span><span class="comment">                          60              150          NS =  10(+)
</span><span class="comment">*</span><span class="comment">                         150              590          NS =  **
</span><span class="comment">*</span><span class="comment">                         590             3000          NS =  64
</span><span class="comment">*</span><span class="comment">                        3000             6000          NS = 128
</span><span class="comment">*</span><span class="comment">                        6000             infinity      NS = 256
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                  (+)  By default some or all matrices of this order 
</span><span class="comment">*</span><span class="comment">                       are passed to the implicit double shift routine
</span><span class="comment">*</span><span class="comment">                       <a name="SLAHQR.199"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> and NS is ignored.  See ISPEC=1 above 
</span><span class="comment">*</span><span class="comment">                       and comments in IPARM for details.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       The asterisks (**) indicate an ad-hoc
</span><span class="comment">*</span><span class="comment">                       function of N increasing from 10 to 64.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            ISPEC=5:  Select structured matrix multiply.
</span><span class="comment">*</span><span class="comment">                      (See <a name="ILAENV.206"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> for details.) Default: 3.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ================================================================
</span><span class="comment">*</span><span class="comment">     Based on contributions by
</span><span class="comment">*</span><span class="comment">        Karen Braman and Ralph Byers, Department of Mathematics,
</span><span class="comment">*</span><span class="comment">        University of Kansas, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ================================================================
</span><span class="comment">*</span><span class="comment">     References:
</span><span class="comment">*</span><span class="comment">       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
</span><span class="comment">*</span><span class="comment">       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
</span><span class="comment">*</span><span class="comment">       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
</span><span class="comment">*</span><span class="comment">       929--947, 2002.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
</span><span class="comment">*</span><span class="comment">       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
</span><span class="comment">*</span><span class="comment">       of Matrix Analysis, volume 23, pages 948--973, 2002.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ================================================================
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Matrices of order NTINY or smaller must be processed by
</span><span class="comment">*</span><span class="comment">     .    <a name="SLAHQR.228"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> because of insufficient subdiagonal scratch space.
</span><span class="comment">*</span><span class="comment">     .    (This is a hard limit.) ====
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== NL allocates some local workspace to help small matrices
</span><span class="comment">*</span><span class="comment">     .    through a rare <a name="SLAHQR.232"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> failure.  NL .GT. NTINY = 11 is
</span><span class="comment">*</span><span class="comment">     .    required and NL .LE. NMIN = <a name="ILAENV.233"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>(ISPEC=1,...) is recom-
</span><span class="comment">*</span><span class="comment">     .    mended.  (The default value of NMIN is 75.)  Using NL = 49
</span><span class="comment">*</span><span class="comment">     .    allows up to six simultaneous shifts and a 16-by-16
</span><span class="comment">*</span><span class="comment">     .    deflation window.  ====
</span><span class="comment">*</span><span class="comment">
</span>      INTEGER            NTINY
      PARAMETER          ( NTINY = 11 )
      INTEGER            NL
      PARAMETER          ( NL = 49 )
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0e0, ONE = 1.0e0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      REAL               HL( NL, NL ), WORKL( NL )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, KBOT, NMIN
      LOGICAL            INITZ, LQUERY, WANTT, WANTZ
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            <a name="ILAENV.253"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      LOGICAL            <a name="LSAME.254"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      EXTERNAL           <a name="ILAENV.255"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>, <a name="LSAME.255"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLACPY.258"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAHQR.258"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>, <a name="SLAQR0.258"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>, <a name="SLASET.258"></a><a href="slaset.f.html#SLASET.1">SLASET</a>, <a name="XERBLA.258"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX, MIN, REAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Decode and check the input parameters. ====
</span><span class="comment">*</span><span class="comment">
</span>      WANTT = <a name="LSAME.267"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'S'</span> )
      INITZ = <a name="LSAME.268"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'I'</span> )
      WANTZ = INITZ .OR. <a name="LSAME.269"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'V'</span> )
      WORK( 1 ) = REAL( MAX( 1, N ) )
      LQUERY = LWORK.EQ.-1
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.<a name="LSAME.274"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .AND. .NOT.WANTT ) THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.276"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'N'</span> ) .AND. .NOT.WANTZ ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
         INFO = -4
      ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
         INFO = -5
      ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
         INFO = -11
      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
         INFO = -13
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Quick return in case of invalid argument. ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="XERBLA.296"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SHSEQR.296"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>'</span>, -INFO )
         RETURN
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( N.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Quick return in case N = 0; nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span>         RETURN
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( LQUERY ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Quick return in case of a workspace query ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLAQR0.309"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
     $                IHI, Z, LDZ, WORK, LWORK, INFO )
<span class="comment">*</span><span class="comment">        ==== Ensure reported workspace size is backward-compatible with
</span><span class="comment">*</span><span class="comment">        .    previous LAPACK versions. ====
</span>         WORK( 1 ) = MAX( REAL( MAX( 1, N ) ), WORK( 1 ) )
         RETURN
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== copy eigenvalues isolated by <a name="SGEBAL.318"></a><a href="sgebal.f.html#SGEBAL.1">SGEBAL</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, ILO - 1
            WR( I ) = H( I, I )
            WI( I ) = ZERO
   10    CONTINUE
         DO 20 I = IHI + 1, N
            WR( I ) = H( I, I )
            WI( I ) = ZERO
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Initialize Z, if requested ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( INITZ )
     $      CALL <a name="SLASET.332"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Z, LDZ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Quick return if possible ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILO.EQ.IHI ) THEN
            WR( ILO ) = H( ILO, ILO )
            WI( ILO ) = ZERO
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== <a name="SLAHQR.342"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>/<a name="SLAQR0.342"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a> crossover point ====
</span><span class="comment">*</span><span class="comment">
</span>         NMIN = <a name="ILAENV.344"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="SHSEQR.344"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>'</span>, JOB( : 1 ) // COMPZ( : 1 ), N, ILO,
     $          IHI, LWORK )
         NMIN = MAX( NTINY, NMIN )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== <a name="SLAQR0.348"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a> for big matrices; <a name="SLAHQR.348"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> for small ones ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( N.GT.NMIN ) THEN
            CALL <a name="SLAQR0.351"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
     $                   IHI, Z, LDZ, WORK, LWORK, INFO )
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Small matrix ====
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLAHQR.357"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
     $                   IHI, Z, LDZ, INFO )
<span class="comment">*</span><span class="comment">
</span>            IF( INFO.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== A rare <a name="SLAHQR.362"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> failure!  <a name="SLAQR0.362"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a> sometimes succeeds
</span><span class="comment">*</span><span class="comment">              .    when <a name="SLAHQR.363"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> fails. ====
</span><span class="comment">*</span><span class="comment">
</span>               KBOT = INFO
<span class="comment">*</span><span class="comment">
</span>               IF( N.GE.NL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Larger matrices have enough subdiagonal scratch
</span><span class="comment">*</span><span class="comment">                 .    space to call <a name="SLAQR0.370"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a> directly. ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLAQR0.372"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR,
     $                         WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Tiny matrices don't have enough subdiagonal
</span><span class="comment">*</span><span class="comment">                 .    scratch space to benefit from <a name="SLAQR0.378"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>.  Hence,
</span><span class="comment">*</span><span class="comment">                 .    tiny matrices must be copied into a larger
</span><span class="comment">*</span><span class="comment">                 .    array before calling <a name="SLAQR0.380"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>. ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLACPY.382"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, N, N, H, LDH, HL, NL )
                  HL( N+1, N ) = ZERO
                  CALL <a name="SLASET.384"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
     $                         NL )
                  CALL <a name="SLAQR0.386"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR,
     $                         WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO )
                  IF( WANTT .OR. INFO.NE.0 )
     $               CALL <a name="SLACPY.389"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, N, N, HL, NL, H, LDH )
               END IF
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Clear out the trash, if necessary. ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
     $      CALL <a name="SLASET.397"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'L'</span>, N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Ensure reported workspace size is backward-compatible with
</span><span class="comment">*</span><span class="comment">        .    previous LAPACK versions. ====
</span><span class="comment">*</span><span class="comment">
</span>         WORK( 1 ) = MAX( REAL( MAX( 1, N ) ), WORK( 1 ) )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== End of <a name="SHSEQR.405"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a> ====
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
